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The  complex  cepstrum  technique  was  investigated  for  possible  use  in  removing 
echoes  caused  by  early  reflections  when  low-frequency,  peaked-response  transducers  are 
calibrated  with  transient  signals.  The  calibration  environment  permits  a high  signal-to- 
noise  ratio  and  a free  choice  of  input  waveforms,  both  of  which  are  advantages  in  com- 
plex cepstrum  processing.  The  proper  method  of  computation  with  this  technique  is 
discussed  in  detail.  Cepstrum  filtering  methods  are  discussed.  Synthetic  measurement  -J— •^^^.7’ 
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20.  ABSTRACT  (Continued) 

data  produced  by  a J9  projector  driven  with  a damped  sinusoidal  pulse  and  a pressure- 
release  reflection  are  shown,  and  real  measurement  data  from  a low-frequency  line 
array  driven  by  a single-cycle  sinusoid  with  complex  surface  and  bottom  reflections  are 
also  shown.  In  both  the  synthetic  and  the  real  data,  the  echoes  were  successfully  re- 
moved using  the  complex  cepstrum  technique.  Data  were  digitized  directly  and  stored 
on  digital  magnetic  tape.  Oversampling  was  reduced  by  sifting  to  achieve  the  correct 
sampling  rate  required  by  the  complex  cepstrum  method.  The  results  indicate  that 
echo  removal  performed  by  complex  cepstrum  processing  can  be  accurate  enough  to 
have  potential  usefulness  in  calibration  procedures. 
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COMPLEX  CEPSTRUM  PROCESSING  OF  DIGITIZED  TR^SIENT  CALIBRATION 
DATA  FOR  REMOVAL  OF  ECHOES 


INTRODUCTION 

Calibration  procedures  employing  broadband  tr^sient  ^g^^.^^^^XsRoto/the 

Stations.  Tha  tunetioa  of  the  ^ 

Sken  here  to  include  such  bodies  of  water  as  lakes  or  pools  as  well  as  tank^  d^torte 

?o“nL“Cz.„«o„  • 

method  under  study  here. 

confusing. 

The  basic  steps  necessary  in  the  complex  cepstrum  process  are  given  in  Fig.  1. 
(Complex  cepstrum  is  the  name  given  to  one  function  in  the  timn  of 
* • fv,  Krtf f r.m  nf  Fie  1 We  shall  refer  to  that  function,  and  sometimes  to  the  entire 

report  is  to  discuss  in  detail  the  computation  method  of  this  system.  A s P P 

to^fmonstrate  the  proper  procedure  for  applying  the  process  to  real  data. 


CEPSTRUM  PROGRAM  THEORY 

A computer  program  that  will  accomplish  the  deconvolution  by  "°"'- 

plex  cepstrum  technique  can  be  written  by  performing  the  steps  called  for  in  F g.  . 

Step  1 consists  of  simply  transforming  the  time  series  into  the  frequency  domain  by 
performing  a fast  Fourier  transform  (FFT)  operation. 
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Fig.  1 — Steps  in  the  complex 
cepstrum  process 


Step  2,  taking  the  complex  logarithm,  is  done  by  employing  the  relation 

log  z = log  r + iO 

where  z = jc  + xy.  But  r,  the  magnitude,  equals + y^,  and  0,  the  phase,  equals 

Log  r and  0 become  the  real  and  imaginary  components  of  a complex  num- 
ber representing  the  complex  logarithm  of  z. 


Steps  3 and  4 are  called  “unwrapping  the  phase.”  They  are  extremely  critical  steps 
in  the  process.  In  order  to  be  able  to  identify  the  p>eriodic  component  representing  the 
echo  in  this  sequence,  it  is  necessary  to  remove  the  numerous  discontinuities  in  the  phase 
that  arise  from  the  multivalued  definition  of  tan’l . Each  time  the  phase  value  reaches 
±7r,  a discontinuity  of  27r  appears  in  the  function.  These  discontinuities  must  be  forced 
out  of  the  function  by  adding  or  subtracting  2ir  to  all  subsequent  values,  as  appropriate. 

This  process  is  performed  in  the  computer  program  be  comparing  each  successive  point 
to  its  neighbor.  If  a difference  of  3 or  more  is  found  between  neighboring  values,  the  se- 
cond point  is  identified  as  a positive  or  negative  discontinuity,  as  the  caise  may  be,  and  an 
adjustment  is  made  to  this  and  all  succeeding  points.  The  value  3 was  arrived  at  empirically. 
The  sampling  used  causes  the  actual  values  obtained  between  discontinuous  points  to  vary. 

If  the  frequency  intervals  are  too  coarse,  discontinuities  can  even  be  lost.  This  can  be 
avoided  by  sampling  for  a long  enough  time  to  provide  adequate  frequency  resolution. 

If  sufficient  time  or  record  length  is  not  available,  zeroes  can  be  added  to  the  time  se- 
quence as  required.  Because  the  phase  data  are  always  antisymmetric,  we  can  generate  the 
second  half  of  the  data  by  forming  an  inverted  mirror  image  of  the  first  half.  Then  we 
need  only  to  phase-unwrap  the  first  half. 
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Removal  of  the  discontinuities  introduces  a huge  linear  slope,  which  is  removed  next, 
in  step  4. 

Improper  or  incomplete  removal  of  the  linear  phase  will  leave  large  components  in 
the  cepstrum  sequence,  which  mask  the  true  amplitude  of  the  echo  impulse  train  com- 
ponents. The  method  we  use  to  remove  the  linear  phase  may  be  explained  by  letting  5 T 
be  the  sampling  interval  and  N be  the  total  number  of  samples.  We  wish  to  eluninate  any 
net  time  delay.  If  there  is  a delay  of  Af  samples,  we  wish  to  shift  the  data  back  by  Mb  T. 
A time  delay  of  MS  T shows  up  in  the  spectrum  as  a multiplicative  factor  of  exp  (iAf5  Tw) 
In  the  phase  of  the  spectrum  (the  imaginary  part  of  the  logarithm)  it  shows  up  as  an  addi- 
tive term  where  = Af6  Tco. 


Now 

u)  = n5  CO 

n = 0,  1,  . . 

and 

2 CO 

_ w fnax 

N 

where 

_ jj 

^ max  ~ 5 T ■ 

So 

2iTMn 

<t>d  - 

n = 0,  1,  . . 

In  this  case,  n ranges  only  to  (lV/2  - 1)  because  the  FFT  requires  that  N be  an  even 
number.  The  phase  function  will  have  a point  of  symmetry  at  N/2,  which  has  a value  of 
zero.  The  total  phase  is  then 

0(«)  = n = 0,  1 ^ 


By  the  integral  definition  of  M,  its  value  can  be  chosen  such  that  the  absolute  value  of 
4>q{NI2)  is  always  less  than  it.  If  Af  is  large,  <»(lV/2)  is  approximately  equal  to  <t)Q{NI2  - 1). 
Then  the  total  phase. 


/ N 

\ /A 

r\  /N 

A . /Af 
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and 


Again,  if  N is  very  large  compared  to  M,  the  last  term  is  negligible.  We  have  seen  that  the 
second  term  must  be  less  than  1.  So  Af  is  the  integer  closest  to[(N/2)-l].  To  form  a 
rounded  value,  set  Af  equal  to  the  integer  value 


of 


and 


for  positive  Af 


for  negative  Af. 


Since  0o(n)  = 0 (n)  - ^^(n),  the  corrected  phase  is  given  by 


•Af> 


00  («)  = 0 (n)  - 'J- 


N 

2/ 


(The  term  in  brackets  is  the  correction  increment,  FAZINC,  used  in  the  computer  program 
given  in  the  Appendix.) 

Step  5,  which  follows,  consists  of  performing  the  inverse  transform  by  an  FFT  oper- 
ation. This  yields  the  complex  cepstrum  (which  is  actually  a real  function  with  no  imag- 
inary part,  because  the  real  part  of  the  log  spectrum  is  symmetric  and  the  imaginary  part 
is  antisymmetric). 

Step  6,  called  linear  filter,  may  apply  to  any  operation  that  is  intended  to  separate 
the  echo  impulse  train  from  the  underlying  cepstrum  curve.  In  our  case  we  wish  to  re- 
move the  echo.  The  simplest  method  might  be  to  zero  all  values  beyond  a chosen  point. 
That  oi)eration  is  called  an  early-pass  filter,  corresponding  to  the  low-pass  filter  in  the  fre- 
quency domain.  Because  the  cepstrum  tends  to  concentrate  the  direct-arrival  signal  in  the 
earliest  part  of  the  curve,  most  of  that  signal  would  be  undisturbed  by  such  an  operation. 
The  process’s  effectiveness  depends  mostly  on  how  well  the  direct-arrival  and  echo  por- 
tions of  the  cepstrum  can  be  separated.  This  in  turn  depends  on  how  close  in  time  the 
echo  and  direct  signal  are. 

A more  specific  method  of  removing  the  echo  would  be  to  identify  the  individual 
points  in  the  echo  impulse  train  and  set  those  to  zero,  leaving  the  intervening  parts  of  the 


4 


NRL  REPORT  8143 


direct-arrival  cepstrum  undisturbed.  This  operation  is  termed  a comb  filter. 
it  allows  remov^  of  echo  signals  of  shorter  delay  than  the  early-pass  filter.  But  zetomg 
i SqueLe  of  points  can  produce  a disturbance  when  the  underlying  cepstrum  is  non  zero 
at  those  points^  A better  method  would  be  to  try  to  replace  the  sequence  of  points  with 
values  from  an  estimate  of  the  echo-free  cepstrum.  In  our  observations  of 
properly  sampled,  band-limited  signals,  the  direct-amval  signal  component  seems  to  have 
L Alternating  point  structure  which  can  be  immediately  visualized  by  refemng  to  it  as 
“Jaws  ” The  best  linear  approximation  for  this  type  of  structure  cannot  ® 
adjacent  points,  but  it  can  be  made  directly  from  the  points  two 

the  central  point.  Thus,  if  an  echo  impulse  point  to  be  removed  is  located  at  point  K 
the  straight-line  approximation  between  K - 2 and  K + 2 is  used  to  ^ f 

ithm  also  works  well  for  smooth  curves.  A general  expression  to  designate  values  of  K is 


X = Af  X « + 1 


„ . 1,  2,  3, 


Where  M is  the  delay  interval  between  signal  and  echo  in  sample  points  and  M equals  the 
total  number  of  sample  points  in  the  sequence.  A second  pass  may  be  formed  by 


X = M X n + 1 


In  both  cases  n is  always  an  integer. 

Steps  7 through  10  are  performed  to  reconstruct  the  time-series  waveform  with  the 
echo  removed.  Step  7 transforms  the  sequence  back  into  the  frequency  domain  with  an 
FFT  operation.  Next,  in  step  8,  the  linear  phase  component  is  restored  using  tl^  same  in- 
crement that  was  determined  in  step  4.  In  step  9 complex  exponentiation  is  performed. 
This  returns  the  spectrum  to  a linear  state. 

Between  steps  1 and  2,  the  first  term  in  the  real  part  of  the  spectrum  (the  DC  tem) 
was  tested  for  algebraic  sign.  This  information  was  stored,  and  if  that  sign  was  found  to 
be  negative  the  sign  of  the  entire  spectrum  was  reversed.  The  original  sign  is  novv  re- 
stored^ just  before  step  10.  Schafer  [21  gives  a good  discussion  of  the  reason  for  this  oper- 
ation. 

Step  10  is  an  inverse  FFT  operation  to  return  the  sequence  to  the  time  domain.  We 
now  have  the  reconstructed  time-series  signal  with  echo  removed. 

With  a computer  program  to  perform  these  steps  in  hand,  we  tried  removing  echoes 
from  some  experimental  data.  The  results  were  generally  unacceptable.  In  order  to  dis- 
cover where  our  problems  were,  we  decided  to  apply  the  program  to  some  more  easUy 
controlled  synthetic  data. 
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SAMPLING  CONSIDERATIONS 

The  first  data  sequence  used  to  test  the  cepstrum  program  was  a simulated  waveform 
generated  from  an  analytic  function  and  forced  to  appear  much  like  the  waveform  used  by 
Ulrych  (Fig.  2).  The  data  sequence  is  128  samples  long.  An  echo  of  0.9  magnitude  and 
delayed  by  12  samples  was  then  added,  making  the  composite  signal  appear  as  in  Fig.  3. 

The  total  wavelet  occupies  about  45  samples.  After  the  first  FFT  (step  2),  one  thing 
immediately  apparent  about  the  simulation  was  the  very  narrow  band-width  of  the  spec- 
trum. This  result  may  be  seen  in  Fig.  4,  where  the  right  half  represents  the  negative-fre- 
quency part  of  the  transform  translated  upwards  (this  is  the  appearance  of  normal  output 
from  an  FFT  routine).  Because  the  energy  in  the  high-frequency  part  of  the  spectrum 
approaches  zero,  the  logarithm  of  the  magnitude  produces  large  negative  values,  and  the 
phase  curve  becomes  erratic.  After  considerable  experimenting,  we  were  able  to  keep  these 
erroneous  values  from  swamping  the  higher  energy  parts  of  the  spectrum  by  introducing  a 
small  impulse,  or  “glitch,”  into  the  time  series.  Its  echo  was  produced  with  the  same  am- 
plitude and  time  delay  as  that  of  the  data  signal,  the  only  intention  being  to  supply  some 
energy  to  the  high-frequency  parts  of  the  spectrum.  When  the  glitch  is  introduced,  the 
real  part  of  the  complex  logeirithm  appears  as  in  Fig.  5.  Note  that  the  effects  of  the  glitch 
and  its  echo  are  apparent  in  the  high-frequency  area.  The  results  of  phase-unwrapping  are 
seen  in  Figs.  6 and  7.  Figure  6 shows  the  imaginary  part  of  the  data.  The  discontinuities 
have  been  removed  from  the  left  half  but  not  from  the  right.  It  was  necessary  to  place  the 
glitch  at  the  center  of  the  data  impulse  to  avoid  introducing  into  the  imaginary  part  another 
phase  component  representing  delay  between  the  data  impulse  and  the  glitch.  This  phase 
component  would  appear  as  a different  slope  over  the  glitch-dominated  frequencies  in  Fig. 

6,  producing  a kink  in  the  linear  ramp.  This  double  slope  would  prevent  the  linear  phase 
from  being  minimized  and  cause  masking  of  the  impulse  train  in  the  cepstrum.  The  com- 
plete unwrapped  phase  signal  appears  in  Fig.  7.  There  is  a considerable  difference  in  ver- 
tical scale  from  Fig.  6.  The  complex  cepstrum  is  seen  in  Fig.  8. 

The  echo  energy  is  represented  by  the  alternating  impulse  train.  The  first  spike  in 
Fig.  8 is  at  the  13th  sample,  and  succeeding  spikes  appear  at  each  12th  sample  following. 
The  right  half  of  Fig.  8 represents  the  same  negative-frequency  translation  effect  shown  in 
Fig.  4,  except  that  the  diminishing  impulse  train  does  not  belong  in  the  negative-frequency 
region.  Furthermore,  if  the  vertical  scale  is  changed  to  magnify  the  curve,  it  can  be  seen 
that  the  impulse  train  wraps  around  the  cepstrum  modulo  128. 

The  linear  filter  operation  was  performed  using  the  special  comb  filter.  The  result  of 
two  passes  of  the  filter,  replacing  every  12th  point  starting  with  13  and  then  every  12th 
point  starting  with  5,  is  shown  in  Fig.  9.  The  result  of  the  reconstructed  time-series  wave- 
form is  shown  in  Fig.  10. 

We  found  no  way  to  obtain  the  preceding  result  with  this  signal  without  resorting  to 
the  glitch  while  still  maintaining  the  number  of  samples  used.  However,  the  problem  posed 
by  near-zero  energy  in  the  high-frequency  part  of  the  spectrum  shown  in  Fig.  4 may  be 
dealt  with  in  a more  satisfactory  way.  By  sampling  at  or  slightly  below  the  Nyquist  rate, 
we  caused  the  important  energy  region  to  just  fill  the  spectrum.  This  results  in  a very  nice- 
looking  cepstrum  which,  unfortunately,  has  an  almost  unrecognizable  time-series  waveform. 
The  method  is  workable,  however,  and  is  the  method  used  in  the  remaining  examples. 
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Fig.  2 — Simulated  waveform  (128  samples) 


Simulated  waveform  plus  echo 
Echo  amplitude  = 0.9 
delay  = 12  samples 


Fig.  4 — Spectrum  of  simulated 
waveform  of  Fig.  3 


Fig.  5 — Real  part  complex  logarithm  of 
simulated  waveform  plus  “glitch” 


Fig.  6 — Imaginary  part  corresponding  to  Pig.  6. 
Left  half:  after  removal  of  phase  discontinuities, 
right  half:  before  removal  of  pha.se  discontinuities. 


Fig.  7 — “Phase  unwrapped”  imaginary  part 
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Fig.  8 — Complex  cepstrum  of  simulated  waveform 


r ig.  9 — Filtered  cepstrum 


Fig.  10  — Reconstructed  simulated 
waveform  minus  echo 


Moreover,  as  will  be  discussed  later,  an  oversampled  waveform  can  always  be  generated 
from  a properly  sampled  one. 

Having  successfully,  if  somewhat  dcvir'usly,  replicated  the  waveforms  of  Ulrych,  we 
began  to  search  for  a better  choice  of  driving  signal.  The  problem  was  obviously  with  the 
sampling  rate;  and  the  solution  was  not  to  put  a glitch  in  the  signal,  but  to  sample  at  ex- 
actly the  right  rate.  One  trial  signal  that  seemed  promising  was  the  impulse  response  of 
an  eight-pole  Chebyshev  filter.  The  idea  here  was  to  have  a sharp  high-frequency  cutoff 
to  avoid  aliasing  problems  when  the  signal  was  sampled  at  exactly  the  correct  rate.  Fig- 
ure 11  shows  the  Chebyshev  filter  impulse  response  as  it  would  appear  in  an  esthetically 
pleasing  but  oversampled  form.  In  Fig.  12  we  see  the  log  spectrum  derived  from  the 
oversampled  signal.  Note  the  sharp  high-frequency  cutoff  and  the  undesirable  low-ampli- 
tude, high-frequency  components.  Figure  13  shows  the  same  waveform  when  the  sampling 
rate  has  been  changed  to  the  proper  value  for  cepstrum  processing;  Fig.  14  is  the  log  spec- 
trum, now  showing  only  about  7-dB  amplitude  veiriation  instead  of  60  dB.  Hardly  any 
aliasing  can  be  seen.  In  Fig.  15  an  echo  has  been  added  to  the  waveform.  This  results  in 
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a cepstrum  (Fig.  16)  which  contains  a very  prominent  sequence  of  single  points  repre- 
senting the  echo.  This  sequence  of  points  was  removed  by  using  the  special  comb  filter, 
and  the  reconstructed  signal  is  shown  in  Fig.  17.  A careful  comparison  of  Fig.  17  and 
Fig.  13  shows  that  the  two  are  virtually  identical. 

When  employing  artificially  generated  signals  as  input  sequences  in  this  manner,  there 
is  a normal  tendency  to  form  the  time-series  waveform  by  adding  to  the  original  signal  a 
reduced-amplitude  replica  delayed  by  some  integral  number  of  sample  points.  In  actual 
practice,  the  echo  will  seldom  be  delayed  by  an  exact  number  of  samples.  The  peak  of 
the  echo  impulse  might  just  as  well  fall  midway  between  two  sample  points.  In  order  to 
investigate  the  effect  of  this  placement  of  echoes  on  the  cepstrum,  we  included  in  the 
original  synthesizing  program  a provision  for  making  the  echo  delay  continuously  variable 
between  sample  points.  Figure  18  shows  the  cepstrum  that  results  when  the  echo  is  intro- 
duced midway  between  sample  points.  Competring  this  figure  with  the  integral  sample  de- 
lay version  (Fig.  16),  we  see  that  the  impulse  train  peaks  have  been  reduced  in  amplitude 
and  spread  so  that  they  are  no  longer  single  points.  .An  attempt  to  remove  the  echo  re- 
presented in  Fig.  18  by  the  comb  filter  is  generally  unsuccessful,  even  with  repeated  passes, 
as  can  be  seen  from  Fig.  19. 

.At  this  point  it  is  clear  that  there  are  two  rules  pertaining  to  sampling  that  must  be 
observed.  First,  the  signal  must  be  sampled  at  exactly  the  proper  rate:  a rate  low  enough 
that  the  spectrum  does  not  extend  beyond  the  frequencies  of  interest,  and  high  enough 
that  no  significant  aliasing  takes  place.  Second,  the  samples  must  be  chosen  so  that  the 
echo  falls  on  or  very  nearly  on  a sample  point.  The  technique  we  used  to  accomplish  this 
was  to  sample  at  as  high  a rate  as  possible  consistent  with  the  other  requirements  of  the 
system  (for  example,  storage).  To  obtain  the  desired  sampling  rate,  we  would  sift  the 
data,  saving  every  Pth  point,  where  P represents  the  divisor  necessary  to  give  proper 
sampling,  and  then  shift  the  starting  point  of  the  sequence  until  the  echo  fell  as  close  as 
(xjssible  to  a sample  point.  With  this  technique  we  were  able  to  obtain  good  results  from 
our  records  of  real  experimental  data. 


Kifi.  11  — Impulse  response  of  eight-pole  Eig.  12  Log  spectrum  of  Chebyshev  filter 

Chrbyshcv  filter 
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3 — Properly  sampled  Chebyshev  filter 
impulse  response  without  echo 


Fig.  14  — Log  spectrum  of  properly 
sampled  waveform 


140  00  0 00 


ig.  15  — Waveform  with  echo  added 


Fig.  16  — Complex  cepstrum  of  signal 
plus  echo 


140  00  0 00 


17  — Recon-structed  waveform  with 
echo  removed 


Fig.  18  — Cepstrum  of  signal  with  echo 
delayed  midway  between  samples 
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OCX)  140  00 

Fig.  19  — Reconstructed  waveform,  mid*sampie  echo 


EXPERIMENT 

In  order  to  test  the  program  with  the  real  experimental  data,  we  used  two  types  of 
signals  that  had  been  recorded  previously.  One  type  of  waveform  was  obtained  by  using 
a damped  sinusoid  as  the  driving  current,  and  the  other  type  by  using  a single-cycle  sine 
wave.  The  first  was  produced  by  driving  a J9  transducer  with  a damped  sinusoidal  wave- 
form and  detecting  with  an  H52  hydrophone.  A reflector,  consisting  of  a large,  square 
piece  of  closed-cell  neoprene  cemented  to  a 1 /4-in. -thick  plexiglass  plate,  had  been  posi- 
tioned so  as  to  supply  an  echo  with  adjustable  delay.  The  hydrophone  signal  was  directly 
digitizetl  using  a Nicolet  digital  oscilloscope  and  recorded  on  a digital  tape  unit.  The  dig- 
itizing rate  was  set  at  the  highest  value  that  would  allow  the  entire  signal  waveform  to  be 
stored  within  the  4096-point  capacity  of  the  oscilloscope.  Results  of  two  of  these  exper- 
iments are  shown  in  Figs.  20  through  27.  The  original  impulses  with  no  echoes  are  repre- 
duced  graphically  in  Figs.  20  and  24.  These  waveforms  were  obtained  by  removing  the 
reflector  plate  from  the  water.  Next,  the  reflector  was  lowered  to  contribute  the  echo 
(Figs.  21  and  25).  The  driving  waveform  in  the  first  four  figures  (Figs.  20-23)  was  a 
damped,  ringing  impul.se  from  an  RLC  circuit.  In  the  second  set  of  figures  (Figs.  24-27), 
the  circuit  was  adjusted  to  give  a much  longer  time  constant,  so  that  only  the  first  half- 
cycle is  seen  here.  The  early  part  of  the  impulse  is  due  to  the  transient  response  of  the 
.19  transducer.  Figures  22  and  26  show  the  cepstrum  resulting  from  each  waveform.  Note 
that  in  the  case  of  a pre.ssure-release  echo,  the  echo  appears  in  the  cepstrum  as  a sequence 
of  negative  impulses  only,  instead  of  as  an  alternating  sequence.  The  same  comb  filter 
that  was  used  previously  has  been  applied;  the  reconstructed  results  are  shown  in  Figs.  23 
and  27. 

Another  type  of  data  used  was  a sample  of  actual  transient  calibration  data  from  an 
element  of  a line  array  hydrophone.  The  data  were  recorded  at  our  Leesburg  Facility. 
Rigging  is  usually  done  there  at  a depth  where  the  pressure-release  surface  reflection  will 
arrive  at  approximately  the  same  time  as  the  hard  reflection  from  the  side  walls  of  the 
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spring.  In  the  case  of  omnidirectional  instruments  this  results  in  a minimum  amplitude 
but  somewhat  distorted  echo.  Figures  28  and  29  show  the  received  waveform  in  the  over- 
sampled and  properly  sampled  forms,  respectively.  Bearing  in  mind  the  distorted  form  of 
the  echo,  we  were  not  surprised  to  see  the  clearly  unspectacular  cepstrum  as  it  appears  in 
Fig.  30.  What  is  surprising  is  the  degree  to  which  the  echo  is  eliminated  with  a single  pass 
of  the  comb  filter  (Fig.  31). 


-r-"  I , . , , , , 

000  14000 

Fig.  20  — Experimentai  waveform  from 
J9  transducer 


Fig.  21  — Experimental  waveform  with  echo 


Fig.  22  — Complex  cepstrum  of 
signal  plus  echo 


Fig.  23  — Reconstructed  waveform 
with  echo  removed 
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y'"" ^ 


Fig.  24  — Experimental  waveform  from 
J9  transducer 


Fig.  26  — Complex  cepstrum  of 
signal  plus  echo 


28000 


Fig.  28  — Transient  calibration  signal 
with  echo,  oversampled 


Fig.  25  — Experimental  waveform  with  echo 


0-00  18000 

Fig.  27  — Reconstructed  waveform  with 
echo  removed 


000  7000 

Fig.  29  — Signal  plus  echo, 
properly  sampled 
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Fig.  30  — Complex  cepstrum  of  transient 
calibration  signal 


Fig.  31  — Reconstructed  signal  with 
echo  removed 


In  all  of  the  previous  examples,  the  properly  sampled  input  signals  and  the  recon- 
structed versions  have  a spiky  appearance,  which  makes  them  difficult  to  relate  to  the 
oversampled  data  we  are  accustomed  to  seeing.  It  is  desirable  to  restore  the  smoothed 
appearance  to  the  sampled  signal  by  interpolation.  The  interpolation  formula  used  was 
derived  from  the  expression; 


sin  2nf^(t  - nT) 

2nf,  (t  - nT) 

n = -» 

which  assumes  that  the  signal  jc(t)  is  band-limited  and  the  spectrum  is  zero  for  I ^ I > /■(, 
and  T = 1/  (2  f^^).  If  interpolation  points  are  chosen  midway  between  sample  points,  the 
computation  is  reduced  to 


x(t) 


= I 


x(nT) 


xl(m  + 1/2)7’]  = 


x(nT) 


(-1)'"  * " 
n{m  + 1/2  - n) 


The  process  of  interpolating  midpoints  is  repeated  until  the  entire  array  of  512  points 
is  filled.  Figures  32  and  33  illustrate  the  results  obtained  using  the  modified  program  on 
the  data  of  Figs.  29  and  31,  respectively. 

Up  to  now  we  have  only  briefly  mentioned  other  obvious  methods  of  linear  filtering 
for  the  cepstrum  (e.g.,  “late-pass”  and  “early-pass”  filtering).  These  are  analogous  to 
high-pass  and  low-pass  in  the  frequency  domain.  In  our  case,  a high  signal-to-noise  ratio 
and  the  ability  to  select  an  optimum  driving  signal  make  it  possible  to  keep  the  cepstrum 
“clean”  enough  for  the  echo  impulse  train  to  be  identified  clearly.  When  the  source  sig- 
nal and  the  echo  are  closely  spaced  in  time,  the  signal  portion  of  the  cepstrum  and  the 
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impulse  train  may  not  be  separated  enough  to  avoid  producing  severe  distortion  in  the 
source  signal  when  early -pass  filtering  is  used.  Stated  another  way,  it  appears  that  the 
tc  ^ ri'ter  is  able  to  remove  echoes  that  have  a shorter  delay  time  than  those  the  early-pass 
1,^  - -^ove.  ’Harly-pass  filtering  can  be  useful  in  cases  where  the  echo  does  not  ap- 

"n  as  a sequence  of  single  points,  either  because  of  distortion  or  because 
u.  etween  sample  points,  as  in  Fig.  18.  Late-pass  filtering  is  of  limited  use- 
ful result  of  reconstruction  indicates  the  time  delay,  which  should  be  de- 

tec'' ...  tpstrum  as  well. 


H ^ T T-  r - T r 1 1 -t T 1 1 1 1 ' ' 

0 00  280  00  0 00  280  00 


Fig.  32  — Transient  calibration  data  with  echo.  Fig.  33  — Reconstructed  signal, 

properly  sampled,  smoothed  echo  removed,  smoothed 

Exponential  weighting,  as  discussed  in  Ulrych  [3]  and  elsewhere,  has  not  been  men- 
tioned in  the  above  discussion.  The  cepstrum  method  does  not  work  at  all  if  the  echo 
sequence  is  not  minimum  phase  (see  Schafer  [2] ).  Exponential  weighting  can  be  useful 
in  forcing  the  echo  components  to  be  minimum  phase,  should  there  be  multiple  echoes. 
Too  much  weighting  can  cause  considerable  inaccuracies  in  recovering  the  source  signal, 
however.  The  effect  of  exponential  weighting  is  chiefly  of  use  when  it  is  more  important 
to  detect  the  delay  time  than  to  recover  the  source  signal. 


CONCLUSION 

We  have  demonstrated  that  echo  removal  performed  by  complex  cepstrum  process- 
ing can  be  accurate  enough  to  have  potential  usefulness  in  calibration  procedures. 
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Appendix  A 

COMPUTER  PROGRAMS  AND  SUBROUTINES 

This  section  contains  listings  of  the  computer  program  and  subroutines  pertinent  to 
this  report.  The  source  language  is  Fortran  IV,  and  all  programs  were  implemented  on  a 
Digital  Equipment  PDP-11/45  computer  with  a multi-user  operating  system.  Subroutines 
that  have  not  been  listed  either  are  specific  to  the  operational  details  of  our  system  and, 
as  such,  would  be  of  little  use  to  others,  or  are  ordinary  utility  subroutines  that  probably 
already  exist  in  a prospective  user’s  library.  In  those  cases,  a full  explanation  of  the  func- 
tional nature  of  each  subroutine  is  given. 
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FORTRAN  IV-PLDS  VW-«4  16:43:13  l»-JAH-77  PAGE  1 

PPTEST.  mf  /ROXTR:  BLOCKS/'VR 

C PROGRAM  PPTESTO)  3 JAR  77 

C 

C 

C THIS  VERSION  USES  SINGLE  PRECISION  FOR 

C 

C ***  EVERYTHING  *** 

C 

C OSES  NO  COMPLEX  FUNCTIONS 

C 

C 


C 

C 

eeoi 

0002 

0003 

c 

0004 

c 

c 

0005 

0006 
0007 
6008 

C 

0000 

0010 

0011 

C 

C 

0012 

0013 

C 

C 

0014 

0015 

6016  4 

0017 

0018 

C 

C 

CO  19 
0020 
0021 
'•-022 

C 

C 

C 

0923 

C 

C 

0024 

C025 

0026 

C027  462 

0028  463 

C 


DIMENSION  IHEA0( 3) , IRAY( 200) , JRAY( 200) , FAKIT(  1024) 

BYTE  ILFA(30) ,FILT 
PI=3. 14159265 

THIS  CALL  SUPPRESSES  FLOATING  ZERO  DIV  ERROR# 

CALL  ERHSET(73, .TRUE. . .FALSE. . .FALSE. , .FALSE. ,) 

OPTIONAL  SUBROUTINE  SUBSTITUTES  ANALYTICALLY  GENERATED 
SIMULATED  WAVEFOHMC AFTER  ULRYCH)  FOR  REAL  DATA. 
VRITE(6,462) 

READ(5.463) lULRY 
IF(  lULRY.EQ. 1)N=I28 
IF(  lULRY.EQ.  DGO  TO  4 

SPECIFY  INPUT  FILE* 

WRITE!  6. 20) 

READ! 5. 30) ICNT, ILFA 
ILFA! ICNT+1)=0 

SET  UP  FILENAME  ARRAY  FOR  INPUT  FILE.  SPECIFY  * OF  PTS 
TO  READ  FROM  FILE. 

WRITE!  6, 142) 

READI5. 109)NF 

SPECIFY  NO  OF  DATA  ARRAY  POINTS  !DP  TO  512).  MUST  BE 
A POWER  OF  2. 

WRITEI6. 106) 

READ! 5, 109) N 
IHR=N*2 
SN= FLOAT!  N) 

ON= FLOAT!  IHN) 

CHOOSE  COMB. EARLY  PASS  OR  LATE  PASS  GATING  TO  BE  DONE 
OH  CEPSTRUM 
WRITE!6. 11) 

READ!5.34)FILT 

IF!  lULRY.Ett.  DCALL  UF!  FAXIT) 

IF!  lULRY.EQ.  DGO  TO  3 

CALL  THE  INPUT  SIGNAL  FILE-READING  SUBROUTINE.  THIS  SUBR. 
ALLOWS  VARIABLE  STARTING  POINT  8 SIFTING  FACTOR. 

LOADS  DATA  IN  ARRAY. 

CALL  FILRED! FAXIT. ILFA, NF) 

SELECT  OPTION  TO  BYPASS  PROGRAM.  GOES  TO  DIRECT 
' INTERPOLATION  OF  SIFTED  INPUT  DATA 
WRITE!  6. 230) 

READ! 5, 463) lORIC 

IF!  lORIG.EQ.  DGO  TO  250 

FORMAT! ’8  IF  YOU  WANT  ULRYCH  TYPE  1:  ’) 

FORMAT!  150) 

OUTPUT  A FILE  FOR  PLOTTING  THE  RAW  DATA 
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FORTRAN 

IV-PLUS 

702-04  16:43:15  19-JAn-77  PACE  2 

PPTEST. 

FT1» 

✓RO^TR: BLOCKS/WR 

0029 

3 

CALL  PLTFILIFAKIT.O.R, 1) 

C 

OPTION  TO  USE  EXPONENTIAL  WEICHTINC 

C 

IF  EXP.VEIGHTING  IS  USED, SPECIFY  EXP., IF  NOT, INPUT  1 

•030 

WIITE(6. 119) 

6031 

REA0(5, 112)GAM 

0032 

IFIGAM.EQ. l.E0)GO  TO  21 

0033 

FAM= 1 . EO^GAM 

0034 

DO  1 1= 1 , IHN,2 

0033 

FAM=FAM*GAM 

0036 

FAK1T(  I)=FAKIT(  I)*FAM 

0037 

1 

CONTINUE 

0038 

10 

FORMAT!  2F9.4) 

C 

TRANSFORM  TO  FREO  DOMAIN 

0039 

21 

CALL  FFTDP(FAKIT,N.-1) 

C 

TRAP  TO  INVERT  SPECTRUM  IF  DC  COMPONENT  IS  NEGATIVE. 

0040 

FAKSIG=FAKIT(  D/ABSCFAKIT!  1>) 

0041 

DO  350  1=1, IHN 

0042 

FAKIT(  I)=FAKIT(  I)#FAKSIG 

0043 

350 

CONTINUE 

C 

OUTPUT  A FILE  FOR  PLOTTIHC  THE  SPECTRUM  DATA. 

0044 

CALL  PLTFIL(FAKIT,0,N, 1) 

C 

TAKE  COMPLEX  LOG  OF  SPECTRUM-THE  HARD  WAY. 

3045 

DO  54  1= 1. IHN,2 

0046 

DBAS=SORT(FAKIT(  I)*FAK;  PC  I)+FAKIT(  I+1)*FAKIT(  I+D) 

0047 

DANG=ATAN2(FAKIT( I+l) ,FAKIT( I)) 

0048 

DOLG=ALOG(DBAS) 

C049 

FAKIT(  I)=DOLG 

0050 

FAKIT( I+1)=DANG 

0051 

54 

CONTINUE 

C 

OUTPUT  A FILE  FOR  PLOTTING  THE  LOG  SPECTRUM  DATA. 

0052 

CALL  PLTFIL(FAKIT,0.N, 1) 

0053 

MP=0 

0054 

NP=0 

C 

UNVRAP  THE  PHASE 

C 

LOOK  FOR  DISCONTINUITIES 

C 

AND  COUNT  POS  8 NEG  ONES  SEPARATELY 

C 

START  « IMAC.PART  OF  2ND  PT.  COMPARE  IT  TO  1ST  PT,  < 

C 

TO  MIDPT. , IF  DIFFERENCE  IS  =>3, IT'S  A DISCONTINUITY 

6055 

DO  51  I = 4,N4^2,2 

0056 

FAZDIF=FAK1T(  I)-FAK1T(  1-2) 

0057 

IF(ABS(FAZD1F)-3.E0)51.41,41 

0058 

41 

IF(FAZD1F)42,999.43 

0059 

43 

HP»NP+1 

0060 

IRAY( NP) = I 

0061 

CO  TO  51 

0062 

42 

MP=MP+1 

3063 

JRAY(MP)  = I 

0064 

51 

CONTINUE 

0065 

IRAY(NP+l)=H+4 

0066 

JRAY(MP+l)  = N+4 

0067 

TDELT=O.EO 

0068 

PDELT=e.EO 

0069 

IFINP.EO.eiCO  TO  305 

C 

UNWRAP  THE  POS  DISCS 

0070 

DO  301  1-l.NP 

0071 

PDELT*  PDELT*2 . EO*P I 
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FORTRAN 

IV-PLOS 

V02-04  16:43:15  19-JA11-77  PAGE  3 

PPTEST. 

FTR 

/'RO/TR;  blocks/vr 

0072 

DO  302  J=IRAY(I),IRAY(l>l)-2,2 

0073 

FAKITl  J) =FAKIT( J) -PDELT 

0074 

302 

CONTINUE 

0075 

301 

CONTINUE 

0076 

305 

IF(MP.EO.0)GO  TO  306 

C 

UNWRAP  THE  NEC  DISCS 

0077 

DO  303  1=1, MP 

0078 

TDELT=TDELT+2 . E0*P I 

0079 

DO  304  J=JRAY( 1) , JRAY(  I+l)-2.2 

0080 

FAKITI J ) = F AKI T( J ) +TDELT 

0081 

304 

CONTINUE 

0082 

303 

CONTINUE 

C 

OUTPUT  A FILE  FOR  PLOTTING  OF  UNWRAPPED  PHASE- 1ST  HALF 

C 

OF  DATA  ONLY. 

0083 

306 

CALL  PLTFIL(FAKIT,O.N, 1) 

C 

TAKE  THE  SLOPE  OUT  OF  1ST  HALF! UNWRAPPED)  OF  IHAG.  DATA 

0084 

DBAS=0.5E0 

0083 

DBAS=SIGH( DBAS.FAKITtN) ) 

0086 

FAZDIF=-P I» INT<  FAKIT(  N) /P I+DBAS) 

0087 

F AZ I NC=  F AZD I F/- ( SN/2 . E0  ) 

0088 

DO  55  1=4, N. 2 

0089 

FI  = FLOAT(  1/2-1) 

0090 

FAKIT(  I)=FAKIT( I)+FAZIRC*FI 

0091 

55 

CONTINUE 

C 

SYNTHESIZE  2ND  HALF  OF  IMAG.DATA  FROM  1ST  HALF. 

0092 

FAKITC H+2)=0E0 

0093 

DO  155  1=4, N, 2 

0094 

FAKIT(  IHN+4-I)=-FAKIT(  I) 

0093 

153 

CONTINUE 

C 

OUTPUT  A FILE  FOR  PLOTTING  OF  COMPLETE  DOCTORED-OP 

C 

PHASE  DATA. 

0096 

CALL  PLTFIL(FAKIT,0,N, 1) 

C 

INVERSE  TRANSFORM  TO  CEPSTRUM  DOMAIN 

0097 

CALL  FFTDPIFAKIT.N, 1) 

C 

OUTPUT  A FILE  FOR  PLOTTING  OF  THE  RAW  CEPSTRUM  DATA. 

0098 

CALL  PLTFILIFAKIT.O.N. 1) 

C 

PERFORM  THE  SELECTED  GATING  OPERATION. 

0099 

IFIFILT.EO. ’A’)GO  TO  130 

0100 

IFIFILT.EO. ’B^GO  TO  140 

0101 

IFIFILT.EO. ’C’lGO  TO  150 

0102 

130 

CALL  COHBD< IHN.FAKIT) 

0103 

CO  TO  160 

0104 

140 

CALL  LOWCDC IHN.FAKIT) 

0103 

GO  TO  160 

0106 

150 

CALL  HI YCD( IHN.FAKIT) 

C 

OUTPUT  A FILE  FOR  PLOTTING  OF  THE  MODIFIED  CEPSTRUM 

0107 

160 

CALL  PLTFIL(FAKIT,0,R, 1) 

C 

SCALE  THE  TRANSFORMED  DATA. 

0108 

DO  12  1=1, IHN,2 

0109 

FAKIT(  I)=FAKIT(  I)/SN 

0110 

12 

CONTINUE 

C 

ZERO  THE  IMAC.PART  OF  CEPSTRUM-JOST  TO  BE  SAFE. 

0111 

DO  52  1=2, IHN,2 

0112 

FAKIT(  I)=O.Ee 

0113 

62 

CONTINUE 

C 

TRANSFORM  BACK  TO  LOG  SPECTRUM! FREQUENCY) 

A4 
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FORTRAN 

IV-PLDS 

V02-04  16:43:15  19-J AN-77  PAGE  4 

PPTEST. 

FTN 

/RO/TR: BLOCKS/NR 

eu4 

CALL  FFTDP(FAKIT.N.-l) 

C 

RESTORE  THE  LINEAR  PHASE  COMPONENT- 1ST  HALF  OF  DATA' 

9115 

DO  56  1=4. N. 2 

0116 

FI=FLOAT( 1/2-1) 

0117 

FAKITI  I)=FAK1T( I)-FAZINC*FI 

0118 

56 

CONTINUE 

C 

SYNTHESIZE  THE  2ND  HALF  OF  I MAC. DATA 

01 19 

FAKITC N+2)=©E0 

0120 

DO  166  1=4, N, 2 

0121 

FAKITl  lHN+4-I)=-FAKlT(  I) 

0122 

166 

CONTINUE 

C 

EXPONENTIATE  PHASE! ANTI LOG) -THE  HARD  WAY. 

0123 

DO  53  1=1, IHN.2 

0124 

DOLG=EXP(FAKIT(  I) ) 

0125 

DANG=COS(FAKIT(  I+l) ) 

0126 

DBAS=SIN(FAKIT(  I+D) 

0127 

FAKIT(  I)=D0LC::!DANG 

0128 

FAKITl  I+1)  = D0LG*DBAS 

0129 

53 

CONTINUE 

C 

REPORT  THE  • OF  POS.SNEG. DISCONTINUITIES  REMOVED  IN 

C 

UNWRAPPING 

0130 

WRITE(6, 109)NP.MP 

0131 

DO  360  1=1, IHN 

C 

RE- INVERT  DATA,  IF  INVERTED  ElARLIER. 

0132 

FAKITl  I)=FAKIT( I)*FAKSIG 

0133 

360 

CONTINUE 

C 

TRANS  BACK  TO  TIME  DATA 

0134 

CALL  FFTDFf FAKIT.N,  1) 

0135 

DO  220  1= 1 , IHN, 2 

0136 

FAKITl  I)  = FAinT(  I)/SN 

0137 

220 

CONTINUE 

C 

REMOVE  THE  EXPONENTIAL  WEIGHTING.  IF  ANY.  THIS  GETS 

C 

STICKY  IN  THE  UPPER  END  OF  ARRAY. 

0138 

FAM= 1E0/CAM 

0139 

IFICAM.EQ. l.E0)GO  TO  250 

0140 

DO  2 1= 1 . IHN, 2 

0141 

FAM=FAM*CAM 

0142 

FAKITl I ) = FAKITl  I ) /FAM 

0143 

2 

CONTINUE 

C 

SMOOTH  THE  UNDERSAMPLED  DATA  WITH  INTERPOLATION.  THIS 

C 

OPERATION  PLACES  INTERPOLATED  POINTS  MIDWAY  BETWEEN 

c 

DATA  POINTS.  REPEATS  THIS  OPERATION  UNTIL  MAX  ARRAY1512) 

c 

SIZE  IS  REACHED. 

0144 

259 

NTRP=SQRT1  FLOATl  5 12/N) ) 

0145 

DO  500  NIT=I,NTRP 

9146 

DO  510  MUL=1.N-1 

0147 

FKSOM=0. 

0148 

SYCN=1-I)**MUL/PI 

9149 

ID=-1 

0150 

CUL=MUL+0.5 

0151 

DO  520  1=1, N 

0152 

SYGN=-SYCN 

0153 

ID=ID+2 

0154 

FKBUM=FKSUM+FAKIT1  ID)*STGN/(CDL-I) 

0155 

520 

CONTINUE 

0156 

FAKITl  2»MUL) =FESim 
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FORTRAN  IV-PLDS 
PPTEST. FTR 


6157 

0158 

6159 

6160 
6161 
6162 
6163 


6164 

6165 

6166 

6167 

6168 

6169 

6170 

6171 

6172 

6173 

6174 

6175 

6176 

6177 

6178 

6179 


510 


530 

500 

C 

C 


11 

20 

30 

34 

106 

109 

111 

112 

1 19 

120 
142 
216 
230 
999 


Ve2-64  16:43:15  19- JAR-77  PAGES 

/RO/TR: BLOCKS/im 

CONTINUE 

N=2*N-1 

IF(NIT.EQ.RTRP)  GO  TO  506 
DO  530  J=1.N 

FAKIT( 2*( N- J)  + 1 ) =FAKIT(  N+ 1-J) 

CONTINUE 

CONTINUE  

OUTPUT  A FILE  FOR  PLOTTING  OF  THE  SMOOTHED  RECONSTRUCTED 
DATA 

CALL  PLTFUL(FAKIT,0,N. 1) 

WRITE! 6, 120) 

FORMAT!  ’ 8C0MB!  A) , EARLY!  B) . LATE!  C) , SELECT  F ILTER:  ’ ) 

FORMAT! 'SF I LE  ’) 

FORMAT!  Q.30A1) 

FORMAT! Al) 

FORMAT! '6 INPUT  NO  SAMP  PTS  *) 

FORMAT!  14) 

FORMAT!  F 16. 5) 

FORMAT!  El 0.6) 

FORMAT! ’8EXP.  WEIGHTING  FACT.  ’) 

FORMAT! ’SDONE’) 

FORMAT! ’SNO. FILE  PTS  TO  READ  ’) 

FORMAT!  2X, I3,2X,2!E16.8,2X) ) 

FORMAT!  ’SFOR  ORIG  DATA  PLOT  ONLY.  TYPE  1 ’) 

END 


PROGRAM  SECTIONS 


NAME 

SIZE 

ATTRIBUTES 

6C0DE1 

664444 

1176 

RO. I.CON.LCL 

6PDATA 

666634 

14 

RO.D.CON.LCL 

6IDATA 

666462 

153 

RW.D.CON.LCL 

6VARS 

611646 

2515 

RV.D.CON.LCL 

6TEMPS 

000616 

4 

RW.D.CON.LCL 

TOTAL  SPACE  ALLOCATED  = 617646  3856 

PPTEST,  PPTEST/'-SP=  PPTEST 
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SUBROUTINE  FFTDP( DATA. NW, ISON) 

THIS  SUBROUTINE  TAKEN  FROH  N.H.  BRENNER  ’THREE  FORTRAN 
PROGRAMS  THAT  PERFORM  THE  COOLEY-TDKEY  FOURIER  TRANSFORM* 
MIT  LINCOLN  LAB  TECH  NOTE  1967-2 
REAL  DATA( 2048) , TEMPR, TEMP I , THETA, VR. WI 
N=2*NN 
J=1 

DO  5 1=1, N. 2 
IF( 1-J) 1.2,2 
TEMPR=DATA( J) 

TEMPI  = DATA(  J+1) 

DATA( J)  = DATA(  I) 

DATA( J+1)=DATA(  I+l) 

DATA( I)=TEMPH 
DATA! I+1)=TEMPI 
M=N/2 

IF(J-M)  5,5,4 
J=J-M 
M=M/2 

IF(M-2)  5,3,3 

j=j+ri 

MMAX=2 

IFIHNAX-N)  7,9,9 
ISTEP=2*rn'LAX 
DO  0 H=l,MnAX,2 
WR= FLOAT! ISCH*(  M-l)  ) 

WI  = FLOAT!  MMAX) 

THETA=3. 14159265Ee*WR/WI 
WR= COS!  THETA) 

Wl  = SIN!  THETA) 

DO  8 l=n.N. ISTEP 
J= I+MMAX 

TEMPR=  WR^DATA! J) -WI*DATA!  J+ 1 ) 

TEMP  I = WR^DATA! J+ 1 ) +WI*DATA!  J) 

DATA!  J) =DATA! I ) -TEMPR 
DATA! J+ 1 ) = DATA! 1+ I ) -TEMP  I 
DATA!  I ) = DATA!  I ) +TEMPR 
DATA!  1+ 1 ) = DATA!  1+ 1 ) +TEMP I 

8 CONTINUE 
MMAX= ISTEP 
GO  TO  6 

9 CONTINUE 
RETURN 
END 


M-*  » noon  „„„  nano 


L.  B.  POCHE,  JR. 


SUBROUTIRE  LOWCD( R. HIDATA) 

THIS  SUBROUTINE  PROVIDES  "EARLY-PASS*  CATIRC.  IT  ZEROES 
ALL  POINTS  SYMMETRICALLY  IN  THE  CEPSTRUM  BETWEEN  THE 
CUTOFF  POINT  AND  THE  CENTRAL  POINT  INCLUDING  THE  CUTOFF 
POINT. 

REAL  HIDATA(N) 

WRITE(6, 1) 

READ(S.2) IRLY 

DO  3 I=IRLY*2+1,N-(2*IRLY+1) .2 

HIDATA( I)=OEO 

CONTINUE 

FORTLATI ’SEARI.Y  PASS  CUTOFF  ’) 

FOHMATC 14) 

lUTURN 

END 


SUBROUTINE  H I YCD(  N, HIDATA) 

THIS  SUBROUTINE  PROVIDES  "LATE-PASS"  CATIRC.  IT  ZEROES 
ALL  POINTS  SYMMETRICALLY  BETWEEN  THE  END  POINTS  OF  THE 
CEJ’STRUM  AND  THE  CUTOFF  POINT.  DOES  HOT  ZERO  THE  CUTOFF 
POINT. 

REAL  HIDATA! N) 

WRITE(6, 1) 

READ! 5.2) LATE 
DO  3 I=I,LATE»:2-3,2 
HIDATA!  I)=OEe 
CONTINUE 

DO  4 l=H-!LATE#2-5) .N.2 
HIDATA! I)=0EO 
CONTINUE 

FORMAT!  • SLATE  PASS  CUTOFF  ’ ) 

FORMAT!  14) 

RETURN 

END 
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C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


I 


i 

1 


7 


6 


9 

8 

1 

2 

3 

4 

5 


SUBROUTINE  COnBD( N. H I DATA) 

THIS  SUBROUTINE  GENERATES  A GATING  COMB  FOR 
REIIOVING  THE  S INGI.E-POINT  COMPONENTS  OF  THE  IMPULSE 
TRAIN.  "DELAY”  SHOULD  BE  ENTERED  AS  THE  NO.  OF  THE 
FIRST  POINT  IN  THE  IMPULSE  TR\IN.  "IN'rERVAL"  SHOULD  BE 
ENTERED  AS  THE  NO.  OF  POINTS  BETViEEH  THE  1ST  fl  2ND  8 
SUBSEQUENT  IMPUI.SE  POINTS.  UP  TO  ID  SEQUENCES  OR 
IMPULSE  TRAIN  VALUFjS  MAY  BE  ENTERED.  CATE:D  POINTS  ARE 
REPLACED  BY  THE  ra;AN  VALUES  OF  THE  NEICHBORING 
POINTS  TWO  POINTS  AWAY  AND  ON  EITHER  SIDE  OF  THE  GATED 
POINT.  STARTING-POINT  SYMMETRY  IS  MAINTAINED  AT  THE'.  ENDS 
OF  niE  CEPSTRUM. 

DIMENSION  IDLY(  10) . IDELI  10) 

REAL  H I DATA( 1 ) . FKO . FK 1 . FK3 . FK4 
BYTE  IANS 
K=  I 

WRrrE(6,  I) 

READ(5.2)  II)LY(  K) 

WB1TE(6,3) 

READ(5,2)  IPF.L(  K) 

WRITE(6.4) 

READ(5.3) IANS 

IF( lANS.EQ. ’N’ )GO  TO  6 

K=K+I 

GO  TO  7 

DC  8 L= I . K 

DO  9 J=2»IDLYtL)-l,H-(2*IDLY(L)-l) ,2*IDEL(L) 

FKl  = HinATA(  J-2) 

FK3-H1DATA<  J+2) 

FKO=HIDATA(  J-4) 

FK4=HIDATA(  J+4) 

H I DATA!  J ) » < FKO-  FK4 ) /'2+FK4 

CONTINUE 

CONTI. MUE 

FORMAT!  • 9DELAY  ’ ) 

KOHJiATC  14) 

FoHfWTI '8  INTERVAL  •) 

FORMAT!  • SMOHE?  ! Y OR  N)  ’) 

FOiUlAT!  Al) 

RETURN 

END 


I 

i 


I 

i 
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SUBROUTINE  UF 


This  subroutine  originally  generated  a 128-point  time  series  function  similar  in  ap- 
liearance  to  Ulrych’s  seismic  wavelet.  It  allowed  keyboard  inputs  of  echo  amplitude,  echo 
delay,  and  glitch  amplitude.  The  expression  used  to  generate  the  waveform  was 


UF{t) 


(t  - 66.5)  (t  - 41.5)e 


81 


(t  - 55.8)^ 
270 


1155 


The  glitch  was  added  at  the  49th  point  and  its  echo  at  49  + delay. 


SUBROUTINE  FILRED 

This  is  a file-reading  subroutine  which  is  used  to  load  input  data  from  the  selected 
file.  It  contains  provision  for  keyboard  inputs  which  specify  the  number  of  points  to  be 
discarded  at  the  beginning  of  the  data  sequence  and  the  sifting  factor  to  be  employed  in 
sampling  the  data.  The  data  are  loaded  into  the  array  in  alternate  positions  correspond- 
ing to  the  real  parts  of  complex  numbers. 


SUBROUTINE  PLTFIL 

This  subroutine  writes  data  from  an  array  of  complex  points  as  an  output  file  to  a 
storage  device  in  a format  compatible  with  the  general-purpose  library  graphics  program  in 
use  at  this  laboratory.  It  provides  for  selection  from  the  keyboard  of  real,  imaginary,  or 
magnitude  data  to  be  written  to  the  output  file. 


SUBROUTINE  PLTFUL 

This  subroutine  is  similar  to  PLTFIL  except  that  is  is  designed  to  output  from  a 
512-point  array  of  real  data  only. 


AlO 


